﻿#include "unapVector.hpp"
#ifdef SW_SLAVE
#include "vectorOps.h"
#endif

namespace UNAP {

template <typename T>
Vector<T>::Vector(Communicator *other_comm)
    : length_(0), values_(NULL), dontDel_(false) {
  commcator_ = other_comm;
}
template <typename T>
Vector<T>::Vector(const label length, Communicator *other_comm)
    : length_(length), values_(NULL), commcator_(other_comm), dontDel_(false) {
  if (length_ > 0) {
    values_ = new T[length_];
    this->SET_zero();
  }
}
template <typename T>
Vector<T>::Vector(const T *val, const label &length, Communicator *other_comm,
                  const bool reUse, const int dontDel)
    : length_(length),
      values_(NULL),
      commcator_(other_comm),
      dontDel_(dontDel) {
  if (reUse) {
    values_ = const_cast<T *>(val);
  } else {
    if (length_ > 0) {
      values_ = new T[length_];
      memcpy(values_, val, length_ * sizeof(T));
    }
  }
}
template <typename T>
Vector<T>::Vector(const label len, const T value, Communicator *other_comm)
    : length_(len), values_(NULL), commcator_(other_comm), dontDel_(false) {
  if (length_ > 0) {
    values_ = new T[length_];
    *this = value;
  }
}

template <typename T>
Vector<T>::Vector() : length_(0), values_(NULL), dontDel_(false) {
  commcator_ = NULL;
}
template <typename T>
Vector<T>::Vector(const label length)
    : length_(length),
      values_(NULL),
      commcator_(&COMM::getGlobalComm()),
      dontDel_(false) {
  if (length_ > 0) {
    values_ = new T[length_];
    this->SET_zero();
  }
}
template <typename T>
Vector<T>::Vector(const T *val, const label &length, const bool reUse,
                  const int dontDel)
    : length_(length),
      values_(NULL),
      commcator_(&COMM::getGlobalComm()),
      dontDel_(dontDel) {
  if (reUse) {
    values_ = const_cast<T *>(val);
  } else {
    if (length_ > 0) {
      values_ = new T[length_];
      memcpy(values_, val, length_ * sizeof(T));
    }
  }
}
template <typename T>
Vector<T>::Vector(const label len, const T value)
    : length_(len),
      values_(NULL),
      commcator_(&COMM::getGlobalComm()),
      dontDel_(false) {
  if (length_ > 0) {
    values_ = new T[length_];
    *this = value;
  }
}

template <typename T>
Vector<T>::Vector(const Vector<T> &v)
    : length_(v.length_),
      values_(NULL),
      commcator_(v.commcator_),
      dontDel_(false) {
  if (length_ > 0) {
    values_ = new T[length_];
    memcpy(values_, v.values_, length_ * sizeof(T));
  }
}
template <typename T>
Vector<T>::Vector(const Array<T> &v)
    : length_(v.size()),
      values_(NULL),
      commcator_(&COMM::getGlobalComm()),
      dontDel_(false) {
  if (length_ > 0) {
    values_ = new T[length_];

    std::copy(v.begin(), v.end(), values_);
  }
}

template <typename T>
Vector<T>::~Vector() {
  if (values_ && !dontDel_) {
    delete[] values_;
    values_ = NULL;
  }
}

template <typename T>
Vector<T> &Vector<T>::operator=(const Vector<T> &v) {
  if (&v == this) {
    return *this;
  } else {
    this->commcator_ = v.commcator_;
    if (length_ != v.length_) {
      if (this->values_) {
        delete[] values_;
        values_ = NULL;
      }

      length_ = v.length_;
      values_ = new T[length_];
    }

    if (length_ > 0) memcpy(values_, v.values_, length_ * sizeof(T));

    return *this;
  }
}
template <>
Vector<scalar> &Vector<scalar>::operator=(const Vector<scalar> &v) {
  if (&v == this) {
    return *this;
  } else {
    this->commcator_ = v.commcator_;
    if (length_ != v.length_) {
      if (this->values_) {
        delete[] values_;
        values_ = NULL;
      }

      length_ = v.length_;
      values_ = new scalar[length_];
    }

    if (length_ > 0) {
      IFNOT_SWACC(length_) {
        memcpy(values_, v.values_, length_ * sizeof(scalar));
      }
#ifdef SW_SLAVE
      else {
        MVM_Arrays arrays1;
        init_MVM_Arrays(&arrays1, length_);
        arrays1.A1Ptr = values_;
        arrays1.A2Ptr = v.values_;
        vectorCopy_host(&arrays1);
      }
#endif
    }
    return *this;
  }
}

template <typename T>
Vector<T> &Vector<T>::operator=(const T &a) {
  if (length_ > 0) forAll(i, length_) values_[i] = a;
  return *this;
}
template <>
Vector<scalar> &Vector<scalar>::operator=(const scalar &a) {
  if (length_ > 0) {
    IFNOT_SWACC(length_) { forAll(i, length_) values_[i] = a; }
#ifdef SW_SLAVE
    else {
      MVM_Arrays arrays1;
      init_MVM_Arrays(&arrays1, length_);
      arrays1.A1Ptr = values_;
      arrays1.k1 = a;
      vectorOps_host(&arrays1, &slave_userFunc_fill);
    }
#endif
  }
  return *this;
}

template <typename T>
Vector<T> &Vector<T>::operator+=(const Vector<T> &v) {
#ifdef DEBUG
  if (length_ != v.length_) {
    commcator_->log() << "ERROR in " << __FILE__ << " " << __LINE__
                      << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }
#endif

  T *val1 = values_;
  const T *val2 = v.values_;
  forAll(i, length_) val1[i] += val2[i];

  return *this;
}
template <>
Vector<scalar> &Vector<scalar>::operator+=(const Vector<scalar> &v) {
#ifdef DEBUG
  if (length_ != v.length_) {
    commcator_->log() << "ERROR in " << __FILE__ << " " << __LINE__
                      << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }
#endif

  IFNOT_SWACC(length_) {
    scalar *val1 = values_;
    const scalar *val2 = v.values_;
    forAll(i, length_) val1[i] += val2[i];
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, length_);
    arrays1.A1Ptr = values_;
    arrays1.A2Ptr = v.values_;
    vectorOps_host(&arrays1, &slave_userFunc_aEaPb);
  }
#endif

  return *this;
}

template <typename T>
Vector<T> &Vector<T>::operator-=(const Vector<T> &v) {
#ifdef DEBUG
  if (length_ != v.length_) {
    commcator_->log() << "ERROR in " << __FILE__ << " " << __LINE__
                      << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }
#endif

  T *val1 = values_;
  const T *val2 = v.values_;
  forAll(i, length_) val1[i] -= val2[i];

  return *this;
}
template <>
Vector<scalar> &Vector<scalar>::operator-=(const Vector<scalar> &v) {
#ifdef DEBUG
  if (length_ != v.length_) {
    commcator_->log() << "ERROR in " << __FILE__ << " " << __LINE__
                      << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }
#endif

  IFNOT_SWACC(length_) {
    scalar *val1 = values_;
    const scalar *val2 = v.values_;
    forAll(i, length_) val1[i] -= val2[i];
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, length_);
    arrays1.A1Ptr = values_;
    arrays1.A2Ptr = v.values_;
    vectorOps_host(&arrays1, &slave_userFunc_aEaMib);
  }
#endif

  return *this;
}

template <typename T>
Vector<T> &Vector<T>::operator+=(const T &a) {
  if (length_ > 0) forAll(i, length_) values_[i] += a;
  return *this;
}
template <>
Vector<scalar> &Vector<scalar>::operator+=(const scalar &a) {
  if (length_ > 0) {
    IFNOT_SWACC(length_) { forAll(i, length_) values_[i] += a; }
#ifdef SW_SLAVE
    else {
      MVM_Arrays arrays1;
      init_MVM_Arrays(&arrays1, length_);
      arrays1.A1Ptr = values_;
      arrays1.k1 = a;
      vectorOps_host(&arrays1, &slave_userFunc_shift);
    }
#endif
  }
  return *this;
}

template <typename T>
Vector<T> &Vector<T>::operator-=(const T &a) {
  if (length_ > 0) forAll(i, length_) values_[i] -= a;
  return *this;
}
template <>
Vector<scalar> &Vector<scalar>::operator-=(const scalar &a) {
  if (length_ > 0) {
    IFNOT_SWACC(length_) { forAll(i, length_) values_[i] -= a; }
#ifdef SW_SLAVE
    else {
      MVM_Arrays arrays1;
      init_MVM_Arrays(&arrays1, length_);
      arrays1.A1Ptr = values_;
      arrays1.k1 = -a;
      vectorOps_host(&arrays1, &slave_userFunc_shift);
    }
#endif
  }
  return *this;
}

template <typename T>
Vector<T> &Vector<T>::operator*=(const T &a) {
  if (length_ > 0) forAll(i, length_) values_[i] *= a;

  return *this;
}
template <>
Vector<scalar> &Vector<scalar>::operator*=(const scalar &a) {
  if (length_ > 0) {
    IFNOT_SWACC(length_) { forAll(i, length_) values_[i] *= a; }
#ifdef SW_SLAVE
    else {
      MVM_Arrays arrays1;
      init_MVM_Arrays(&arrays1, length_);
      arrays1.A1Ptr = values_;
      arrays1.k1 = a;
      vectorOps_host(&arrays1, &slave_userFunc_scale);
    }
#endif
  }
  return *this;
}

template <typename T>
Vector<T> &Vector<T>::operator/=(const T &a) {
  if (length_ > 0) forAll(i, length_) values_[i] /= a;
  return *this;
}
template <>
Vector<scalar> &Vector<scalar>::operator/=(const scalar &a) {
  if (length_ > 0) {
    IFNOT_SWACC(length_) { forAll(i, length_) values_[i] /= a; }
#ifdef SW_SLAVE
    else {
      MVM_Arrays arrays1;
      init_MVM_Arrays(&arrays1, length_);
      arrays1.A1Ptr = values_;
      arrays1.k1 = 1. / a;
      vectorOps_host(&arrays1, &slave_userFunc_scale);
    }
#endif
  }
  return *this;
}

template <typename T>
const Vector<T> Vector<T>::operator+(const Vector<T> &v) const {
  return Vector<T>(*this) += v;
}

template <typename T>
const Vector<T> Vector<T>::operator-(const Vector<T> &v) const {
  return Vector<T>(*this) -= v;
}

template <typename T>
const Vector<T> Vector<T>::operator*(const T &a) const {
  return Vector<T>(*this) *= a;
}

template <typename T>
void Vector<T>::SET_size(const label newSize) {
#ifdef DEBUG
  if (newSize < 0) {
    commcator_->log() << "Error in vector SET_size: "
                      << "bad new set size " << newSize << ENDL;

    ERROR_EXIT;
  }
#endif
  if (newSize == 0) {
    this->length_ = 0;
    if (this->values_) {
      delete[] this->values_;
      this->values_ = NULL;
    }
  } else if (newSize != this->length_) {
    T *nv = new T[newSize];

    if (this->length_) {
      label i = MIN(this->length_, newSize);
      T *vv = &this->values_[i];
      T *av = &nv[i];
      while (i--) *--av = *--vv;
    }

    if (this->values_) {
      delete[] this->values_;
      this->values_ = NULL;
    }

    this->length_ = newSize;
    this->values_ = nv;
  }
}

template <typename T>
void Vector<T>::scale(T a) {
  forAll(i, length_) values_[i] *= a;
}
template <>
void Vector<scalar>::scale(scalar a) {
  IFNOT_SWACC(length_) { forAll(i, length_) values_[i] *= a; }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, length_);
    arrays1.A1Ptr = values_;
    arrays1.k1 = a;
    vectorOps_host(&arrays1, &slave_userFunc_scale);
  }
#endif
}

template <typename T>
T Vector<T>::Sum() const {
  T sum = 0;
  forAll(i, length_) sum += values_[i];
  reduceSum(&sum, commcator_);
  return sum;
}
template <>
scalar Vector<scalar>::Sum() const {
  scalar sum = 0;

  IFNOT_SWACC(length_) { forAll(i, length_) sum += values_[i]; }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, length_);
    arrays1.A2Ptr = values_;
    arrays1.k1Ptr = &sum;
    gSum_host(&arrays1, &slave_userFunc_sum);
  }
#endif

  reduceSum(&sum, commcator_);
  return sum;
}

template <typename T>
T Vector<T>::SumMag() const {
  T sum = 0;
  forAll(i, length_) sum += ABS(values_[i]);
  reduceSum(&sum, commcator_);
  return sum;
}
template <>
scalar Vector<scalar>::SumMag() const {
  scalar sum = 0;

  IFNOT_SWACC(length_) { forAll(i, length_) sum += ABS(values_[i]); }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, length_);
    arrays1.A2Ptr = values_;
    arrays1.k1Ptr = &sum;
    gSum_host(&arrays1, &slave_userFunc_sumMag);
  }
#endif

  reduceSum(&sum, commcator_);

  return sum;
}

template <typename T>
T Vector<T>::SumSqr() const {
  T sum = 0;
  forAll(i, length_) sum += values_[i] * values_[i];
  reduceSum(&sum, commcator_);
  return sum;
}
template <>
scalar Vector<scalar>::SumSqr() const {
  scalar sum = 0;

  IFNOT_SWACC(length_) { forAll(i, length_) sum += values_[i] * values_[i]; }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, length_);
    arrays1.A2Ptr = values_;
    arrays1.k1Ptr = &sum;
    gSum_host(&arrays1, &slave_userFunc_sumSqr);
  }
#endif
  reduceSum(&sum, commcator_);
  return sum;
}

template <typename T>
T Vector<T>::SumSqrt() const {
  return sqrt(SumSqr());
}

template <typename T>
void Vector<T>::SET_zero() {
  *this = (T)0;
}

template <typename T>
void reciprocal(Vector<T> &y) {
  const label rowSize = y.size();
  T *yP = y.data();
  forAll(i, rowSize) { yP[i] = 1.0 / yP[i]; }
}
template <>
void reciprocal(Vector<scalar> &y) {
  const label rowSize = y.size();
  IFNOT_SWACC(rowSize) {
    scalar *yP = y.data();
    forAll(i, rowSize) { yP[i] = 1.0 / yP[i]; }
  }
#ifdef SW_SLAVE
  else {
    const Vector<scalar> tmp(y);
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, rowSize);
    arrays1.A1Ptr = y.data();
    arrays1.A2Ptr = tmp.data();
    vectorOps_host(&arrays1, &slave_userFunc_aE1Db);
  }
#endif
}

template <typename T>
void reciprocal(Vector<T> &v1, const Vector<T> &v2) {
#ifdef DEBUG
  if (v1.size() != v2.size()) {
    v1.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label rowSize = v1.size();
  T *v1P = v1.data();
  const T *v2P = v2.data();
  forAll(i, rowSize) { v1P[i] = 1.0 / v2P[i]; }
}
template <>
void reciprocal(Vector<scalar> &v1, const Vector<scalar> &v2) {
#ifdef DEBUG
  if (v1.size() != v2.size()) {
    v1.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label rowSize = v1.size();
  IFNOT_SWACC(rowSize) {
    scalar *v1P = v1.data();
    const scalar *v2P = v2.data();
    forAll(i, rowSize) { v1P[i] = 1.0 / v2P[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, rowSize);
    arrays1.A1Ptr = v1.data();
    arrays1.A2Ptr = v2.data();
    vectorOps_host(&arrays1, &slave_userFunc_aE1Db);
  }
#endif
}

template <typename T>
T dot(const Vector<T> &v1, const Vector<T> &v2) {
#ifdef DEBUG
  if (v1.size() != v2.size()) {
    v1.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  T res = 0;
  const T *v1P = v1.data(), *v2P = v2.data();
  const label size = v1.size();
  forAll(i, size) { res += v1P[i] * v2P[i]; }
  reduceSum(&res, v1.getCommunicator());
  return res;
}
template <>
scalar dot(const Vector<scalar> &v1, const Vector<scalar> &v2) {
#ifdef DEBUG
  if (v1.size() != v2.size()) {
    v1.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  scalar res = 0;
  const label size = v1.size();
  IFNOT_SWACC(size) {
    const scalar *v1P = v1.data(), *v2P = v2.data();
    forAll(i, size) { res += v1P[i] * v2P[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, size);
    arrays1.A2Ptr = v1.data();
    arrays1.A3Ptr = v2.data();
    arrays1.k1Ptr = &res;
    gSum_host(&arrays1, &slave_userFunc_dot);
  }
#endif

  reduceSum(&res, v1.getCommunicator());
  return res;
}

template <typename T>
void axpby(Vector<T> &y, T a, T b, const Vector<T> &x) {
#ifdef DEBUG
  if (x.size() != y.size()) {
    y.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  T *yP = y.data();
  const T *xP = x.data();
  const label size = y.size();
  forAll(i, size) { yP[i] = a * xP[i] + b * yP[i]; }
}
template <>
void axpby(Vector<scalar> &y, scalar a, scalar b, const Vector<scalar> &x) {
#ifdef DEBUG
  if (x.size() != y.size()) {
    y.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label size = y.size();
  IFNOT_SWACC(size) {
    scalar *yP = y.data();
    const scalar *xP = x.data();
    forAll(i, size) { yP[i] = a * xP[i] + b * yP[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, size);
    arrays1.A1Ptr = y.data();
    arrays1.A2Ptr = x.data();
    arrays1.k1 = a;
    arrays1.k2 = b;
    vectorOps_host(&arrays1, &slave_userFunc_axpby);
  }
#endif
}

template <typename T>
void axpy(Vector<T> &y, T a, const Vector<T> &x) {
#ifdef DEBUG
  if (x.size() != y.size()) {
    y.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  T *yP = y.data();
  const T *xP = x.data();
  const label size = y.size();
  forAll(i, size) { yP[i] = a * xP[i] + yP[i]; }
}
template <>
void axpy(Vector<scalar> &y, scalar a, const Vector<scalar> &x) {
#ifdef DEBUG
  if (x.size() != y.size()) {
    y.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label size = y.size();
  IFNOT_SWACC(size) {
    scalar *yP = y.data();
    const scalar *xP = x.data();
    forAll(i, size) { yP[i] = a * xP[i] + yP[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, size);
    arrays1.A1Ptr = y.data();
    arrays1.A2Ptr = x.data();
    arrays1.k1 = a;
    vectorOps_host(&arrays1, &slave_userFunc_axpy);
  }
#endif
}

template <typename T>
void aypx(Vector<T> &y, T a, const Vector<T> &x) {
#ifdef DEBUG
  if (x.size() != y.size()) {
    y.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  T *yP = y.data();
  const T *xP = x.data();
  const label size = y.size();
  forAll(i, size) { yP[i] = xP[i] + a * yP[i]; }
}
template <>
void aypx(Vector<scalar> &y, scalar a, const Vector<scalar> &x) {
#ifdef DEBUG
  if (x.size() != y.size()) {
    y.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of two vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label size = y.size();
  IFNOT_SWACC(size) {
    scalar *yP = y.data();
    const scalar *xP = x.data();
    forAll(i, size) { yP[i] = xP[i] + a * yP[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, size);
    arrays1.A1Ptr = y.data();
    arrays1.A2Ptr = x.data();
    arrays1.k1 = a;
    vectorOps_host(&arrays1, &slave_userFunc_aypx);
  }
#endif
}

template <typename T>
void waxpy(Vector<T> &w, T a, const Vector<T> &x, const Vector<T> &y) {
#ifdef DEBUG
  if (x.size() != y.size() || w.size() != x.size()) {
    w.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of  vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif
  T *wP = w.data();
  const label rowSize = w.size();
  const T *xP = x.data(), *yP = y.data();

  forAll(i, rowSize) { wP[i] = a * xP[i] + yP[i]; }
}
template <>
void waxpy(Vector<scalar> &w, scalar a, const Vector<scalar> &x,
           const Vector<scalar> &y) {
#ifdef DEBUG
  if (x.size() != y.size() || w.size() != x.size()) {
    w.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of  vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif
  const label rowSize = w.size();
  IFNOT_SWACC(rowSize) {
    scalar *wP = w.data();
    const scalar *xP = x.data(), *yP = y.data();
    forAll(i, rowSize) { wP[i] = a * xP[i] + yP[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, rowSize);
    arrays1.A1Ptr = w.data();
    arrays1.A2Ptr = x.data();
    arrays1.A3Ptr = y.data();
    arrays1.k1 = a;
    vectorOps_host(&arrays1, &slave_userFunc_waxpy);
  }
#endif
}

template <typename T>
void axpbypz(Vector<T> &y, T a, T b, const Vector<T> &x, const Vector<T> &z) {
#ifdef DEBUG
  if (x.size() != y.size() || z.size() != x.size()) {
    y.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of  vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  T *yP = y.data();
  const label rowSize = y.size();
  const T *xP = x.data(), *zP = z.data();
  forAll(i, rowSize) { yP[i] = a * xP[i] + b * yP[i] + zP[i]; }
}
template <>
void axpbypz(Vector<scalar> &y, scalar a, scalar b, const Vector<scalar> &x,
             const Vector<scalar> &z) {
#ifdef DEBUG
  if (x.size() != y.size() || z.size() != x.size()) {
    y.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of  vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label rowSize = y.size();
  IFNOT_SWACC(rowSize) {
    scalar *yP = y.data();
    const scalar *xP = x.data(), *zP = z.data();
    forAll(i, rowSize) { yP[i] = a * xP[i] + b * yP[i] + zP[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, rowSize);
    arrays1.A1Ptr = y.data();
    arrays1.A2Ptr = x.data();
    arrays1.A3Ptr = z.data();
    arrays1.k1 = a;
    arrays1.k2 = b;
    vectorOps_host(&arrays1, &slave_userFunc_axpbypz);
  }
#endif
}

template <typename T>
void pointwiseMul(Vector<T> &w, const Vector<T> &x, const Vector<T> &y) {
#ifdef DEBUG
  if (x.size() != y.size() || w.size() != x.size()) {
    w.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of  vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label rowSize = w.size();
  T *wP = w.data();
  const T *xP = x.data();
  const T *yP = y.data();
  forAll(i, rowSize) { wP[i] = xP[i] * yP[i]; }
}
template <>
void pointwiseMul(Vector<scalar> &w, const Vector<scalar> &x,
                  const Vector<scalar> &y) {
#ifdef DEBUG
  if (x.size() != y.size() || w.size() != x.size()) {
    w.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of  vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label rowSize = w.size();
  IFNOT_SWACC(rowSize) {
    scalar *wP = w.data();
    const scalar *xP = x.data();
    const scalar *yP = y.data();
    forAll(i, rowSize) { wP[i] = xP[i] * yP[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, rowSize);
    arrays1.A1Ptr = w.data();
    arrays1.A2Ptr = x.data();
    arrays1.A3Ptr = y.data();

    vectorOps_host(&arrays1, &slave_userFunc_pointwiseMul);
  }
#endif
}

template <typename T>
void pointwiseDiv(Vector<T> &w, const Vector<T> &x, const Vector<T> &y) {
#ifdef DEBUG
  if (x.size() != y.size() || w.size() != x.size()) {
    w.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of  vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label rowSize = w.size();
  T *wP = w.data();
  const T *xP = x.data();
  const T *yP = y.data();
  forAll(i, rowSize) { wP[i] = xP[i] / yP[i]; }
}
template <>
void pointwiseDiv(Vector<scalar> &w, const Vector<scalar> &x,
                  const Vector<scalar> &y) {
#ifdef DEBUG
  if (x.size() != y.size() || w.size() != x.size()) {
    w.getCommunicator()->log()
        << "ERROR in " << __FILE__ << " " << __LINE__
        << ": The length of  vectors is not same!" << ENDL;
    ERROR_EXIT;
  }

#endif

  const label rowSize = w.size();
  IFNOT_SWACC(rowSize) {
    scalar *wP = w.data();
    const scalar *xP = x.data();
    const scalar *yP = y.data();
    forAll(i, rowSize) { wP[i] = xP[i] / yP[i]; }
  }
#ifdef SW_SLAVE
  else {
    MVM_Arrays arrays1;
    init_MVM_Arrays(&arrays1, rowSize);
    arrays1.A1Ptr = w.data();
    arrays1.A2Ptr = x.data();
    arrays1.A3Ptr = y.data();

    vectorOps_host(&arrays1, &slave_userFunc_pointwiseDiv);
  }
#endif
}

template class Vector<scalar>;
template class Vector<label32>;
template class Vector<label64>;
template class Vector<char>;
template class Vector<std::complex<scalar> >;

template void reciprocal(Vector<label> &y);

template void reciprocal(Vector<label> &y, const Vector<label> &x);

template label dot(const Vector<label> &v1, const Vector<label> &v2);

template void axpby(Vector<label> &y, label a, label b, const Vector<label> &x);

template void axpy(Vector<label> &y, label a, const Vector<label> &x);

template void aypx(Vector<label> &y, label a, const Vector<label> &x);

template void waxpy(Vector<label> &w, label a, const Vector<label> &x,
                    const Vector<label> &y);

template void axpbypz(Vector<label> &y, label a, label b,
                      const Vector<label> &x, const Vector<label> &z);

template void pointwiseMul(Vector<label> &w, const Vector<label> &x,
                           const Vector<label> &y);

template void pointwiseDiv(Vector<label> &w, const Vector<label> &x,
                           const Vector<label> &y);

}  // namespace UNAP
